/*
 * #%L
 * Osm2garminAPI
 * %%
 * Copyright (C) 2011 Frantisek Mantlik <frantisek at mantlik.cz>
 * %%
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as
 * published by the Free Software Foundation, either version 2 of the 
 * License, or (at your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public 
 * License along with this program.  If not, see
 * <http://www.gnu.org/licenses/gpl-2.0.html>.
 * #L%
 */
package be.dch.srtm2osm.domain;

//----------------------------------------------------------

import java.util.ArrayList;
import java.util.List;

//----------------------------------------------------------
public class Contours {

    // Below, constant data members:
    // private final static boolean SHOW_NUMBERS = true;
    private double zMaxMin = 1.0E+10;
    // private double zMinMin = -1.0E+10;
    // Below, data members which store the grid steps,
    // the z values, the interpolation flag, the dimensions
    // of the contour plot and the increments in the grid:
    private final int xSteps, ySteps;
    private final float z[][];
    // boolean logInterpolation = false; // not implemented
    private final double xmin, ymin, deltaX, deltaY;
    private List<Contour> contours;
    // Below, data members, most of which are adapted from
    // Fortran variables in Snyder's code:
    private int ncv;
    private double cv[] = new double[2];
    private final double contourStep;
    private final int l1[] = new int[4];
    private final int l2[] = new int[4];
    private final int ij[] = new int[2];
    private final int i1[] = new int[2];
    private final int i2[] = new int[2];
    private final int i3[] = new int[6];
    private int ibkey, icur, jcur, ii, jj, elle, ix, iedge, iflag, ni, ks;
    private int cntrIndex, prevIndex;
    private int idir, nxidir, k;
    private double z1, z2, cval, zMax, zMin;
    private final double intersect[] = new double[4];
    private final double xy[] = new double[2];
    private final double prevXY[] = new double[2];
    private boolean jump;

    // -------------------------------------------------------
    // A constructor method.
    // -------------------------------------------------------
    public Contours(final float[][] z, final int nx, final int ny, final double xmin, final double ymin, final double step, final double contourStep, final double nanValue) {
	this.z = z;
	this.xSteps = nx;
	this.ySteps = ny;
	this.xmin = xmin;
	this.ymin = ymin;
	this.deltaX = step;
	this.deltaY = step;
	this.contourStep = contourStep;
	this.zMaxMin = nanValue - 1;
	// zMinMin = -zMaxMin;
    }

    public List<Contour> makeContours() {
	this.assignContourValues();
	final int workLength = 2 * this.xSteps * this.ySteps * this.ncv;
	boolean workSpace[]; // Allocate below if data valid
	this.contours = new ArrayList<Contour>();
	/* try { */
	if (this.cv[0] != this.cv[1]) { // Valid data
	    workSpace = new boolean[workLength];
	    this.contourPlotKernel(workSpace);
	}
	/*
	 * } catch(Exception ex) { ex.printStackTrace(); return null; }
	 */
	return this.contours;
    }

    // -------------------------------------------------------
    private int sign(final int a, final int b) {
	int c = Math.abs(a);
	if (b < 0) {
	    c = -c;
	}
	return c;
    }

    // -------------------------------------------------------
    // "invalidData" sets the first two components of the
    // contour value array to equal values, thus preventing
    // subsequent drawing of the contour plot.
    // -------------------------------------------------------
    private void invalidData() {
	this.cv[0] = (float) 0.0;
	this.cv[1] = (float) 0.0;
    }

    // -------------------------------------------------------
    // "GetExtremes" scans the data in "z" in order
    // to assign values to "zMin" and "zMax".
    // -------------------------------------------------------
    private void getExtremes() {
	int i, j;
	double here;
	boolean dataSet = false;
	for (i = 0; i < this.xSteps; i++) {
	    for (j = 0; j < this.ySteps; j++) {
		here = this.z[i][j];
		if (here < this.zMaxMin) {
		    if (dataSet) {
			if (this.zMin > here) {
			    this.zMin = here;
			}
			if (this.zMax < here) {
			    this.zMax = here;
			}
		    } else {
			this.zMin = here;
			this.zMax = here;
			dataSet = true;
		    }
		}
	    }
	}
	if (this.zMin == this.zMax) {
	    this.invalidData();
	}
    }

    // -------------------------------------------------------
    // "AssignContourValues" interpolates between "zMin" and
    // "zMax", either logarithmically or linearly, in order
    // to assign contour values to the array "cv".
    // -------------------------------------------------------
    private void assignContourValues() {
	if (this.contourStep <= 0d) {
	    this.invalidData();
	} else {
	    int i;
	    this.getExtremes();
	    /*
	     * if ((logInterpolation) && (zMin <= 0.0)) { invalidData(); } if
	     * (logInterpolation) { double temp = Math.log(zMin);
	     * 
	     * delta = (Math.log(zMax) - temp) / ncv; for (i = 0; i < ncv; i++)
	     * { cv[i] = (float) Math.exp(temp + (i + 1) * delta); } } else {
	     */
	    final double firstContour = (Math.floor(this.zMin / this.contourStep) * this.contourStep) + this.contourStep;
	    this.ncv = (int) Math.floor((this.zMax - firstContour) / this.contourStep) + 1;
	    if (this.ncv <= 1) {
		this.invalidData();
	    } else {
		this.cv = new double[this.ncv];
		for (i = 0; i < this.ncv; i++) {
		    this.cv[i] = firstContour + (i * this.contourStep);
		}
	    }
	}
	// }
    }

    // -------------------------------------------------------
    // "drawKernel" is the guts of drawing and is called
    // directly or indirectly by "contourPlotKernel" in order
    // to draw a segment of a contour or to set the pen
    // position "prevXY". Its action depends on "iflag":
    //
    // iflag == 1 means Continue a contour
    // iflag == 2 means Start a contour at a boundary
    // iflag == 3 means Start a contour not at a boundary
    // iflag == 4 means Finish contour at a boundary
    // iflag == 5 means Finish closed contour not at boundary
    // iflag == 6 means Set pen position
    //
    // If the constant "SHOW_NUMBERS" is true then when
    // completing a contour ("iflag" == 4 or 5) the contour
    // index is drawn adjacent to where the contour ends.
    // -------------------------------------------------------
    private void drawKernel() {
	double prevU, prevV, u, v;

	if ((this.iflag == 1) || (this.iflag == 4) || (this.iflag == 5)) {
	    if (this.cntrIndex != this.prevIndex) { // Must change colour
		// newContour();
		this.prevIndex = this.cntrIndex;
	    }
	    prevU = ((this.prevXY[0] - 1.0) * this.deltaX) + this.xmin;
	    prevV = ((this.prevXY[1] - 1.0) * this.deltaY) + this.ymin;
	    u = ((this.xy[0] - 1.0) * this.deltaX) + this.xmin;
	    v = ((this.xy[1] - 1.0) * this.deltaY) + this.ymin;

	    if ((this.z1 < this.zMaxMin) && (this.z2 < this.zMaxMin)) {
		this.contourSegment(prevU, prevV, u, v, this.cv[this.cntrIndex]);
	    }
	}
	this.prevXY[0] = this.xy[0];
	this.prevXY[1] = this.xy[1];
    }

    // -------------------------------------------------------
    // "detectBoundary"
    // -------------------------------------------------------
    private void detectBoundary() {
	this.ix = 1;
	if (this.ij[1 - this.elle] != 1) {
	    this.ii = this.ij[0] - this.i1[1 - this.elle];
	    this.jj = this.ij[1] - this.i1[this.elle];
	    if (this.z[this.ii - 1][this.jj - 1] <= this.zMaxMin) {
		this.ii = this.ij[0] + this.i2[this.elle];
		this.jj = this.ij[1] + this.i2[1 - this.elle];
		if (this.z[this.ii - 1][this.jj - 1] < this.zMaxMin) {
		    this.ix = 0;
		}
	    }
	    if (this.ij[1 - this.elle] >= this.l1[1 - this.elle]) {
		this.ix = this.ix + 2;
		return;
	    }
	}
	this.ii = this.ij[0] + this.i1[1 - this.elle];
	this.jj = this.ij[1] + this.i1[this.elle];
	if (this.z[this.ii - 1][this.jj - 1] > this.zMaxMin) {
	    this.ix = this.ix + 2;
	    return;
	}
	if (this.z[this.ij[0]][this.ij[1]] >= this.zMaxMin) {
	    this.ix = this.ix + 2;
	}
    }

    // -------------------------------------------------------
    // "routineLabel020" corresponds to a block of code
    // starting at label 20 in Synder's subroutine "GCONTR".
    // -------------------------------------------------------
    private boolean routineLabel020() {
	this.l2[0] = this.ij[0];
	this.l2[1] = this.ij[1];
	this.l2[2] = -this.ij[0];
	this.l2[3] = -this.ij[1];
	this.idir = 0;
	this.nxidir = 1;
	this.k = 1;
	this.ij[0] = Math.abs(this.ij[0]);
	this.ij[1] = Math.abs(this.ij[1]);
	if (this.z[this.ij[0] - 1][this.ij[1] - 1] > this.zMaxMin) {
	    this.elle = this.idir % 2;
	    this.ij[this.elle] = this.sign(this.ij[this.elle], this.l1[this.k - 1]);
	    return true;
	}
	this.elle = 0;
	return false;
    }

    // -------------------------------------------------------
    // "routineLabel050" corresponds to a block of code
    // starting at label 50 in Synder's subroutine "GCONTR".
    // -------------------------------------------------------
    private boolean routineLabel050() {
	while (true) {
	    if (this.ij[this.elle] >= this.l1[this.elle]) {
		if (++this.elle <= 1) {
		    continue;
		}
		this.elle = this.idir % 2;
		this.ij[this.elle] = this.sign(this.ij[this.elle], this.l1[this.k - 1]);
		if (this.routineLabel150()) {
		    return true;
		}
		continue;
	    }
	    this.ii = this.ij[0] + this.i1[this.elle];
	    this.jj = this.ij[1] + this.i1[1 - this.elle];
	    if (this.z[this.ii - 1][this.jj - 1] > this.zMaxMin) {
		if (++this.elle <= 1) {
		    continue;
		}
		this.elle = this.idir % 2;
		this.ij[this.elle] = this.sign(this.ij[this.elle], this.l1[this.k - 1]);
		if (this.routineLabel150()) {
		    return true;
		}
		continue;
	    }
	    break;
	}
	this.jump = false;
	return false;
    }

    // -------------------------------------------------------
    // "routineLabel150" corresponds to a block of code
    // starting at label 150 in Synder's subroutine "GCONTR".
    // -------------------------------------------------------
    private boolean routineLabel150() {
	while (true) {
	    // ------------------------------------------------
	    // Lines from z[ij[0]-1][ij[1]-1]
	    // to z[ij[0] ][ij[1]-1]
	    // and z[ij[0]-1][ij[1]]
	    // are not satisfactory. Continue the spiral.
	    // ------------------------------------------------
	    if (this.ij[this.elle] < this.l1[this.k - 1]) {
		this.ij[this.elle]++;
		if (this.ij[this.elle] > this.l2[this.k - 1]) {
		    this.l2[this.k - 1] = this.ij[this.elle];
		    this.idir = this.nxidir;
		    this.nxidir = this.idir + 1;
		    this.k = this.nxidir;
		    if (this.nxidir > 3) {
			this.nxidir = 0;
		    }
		}
		this.ij[0] = Math.abs(this.ij[0]);
		this.ij[1] = Math.abs(this.ij[1]);
		if (this.z[this.ij[0] - 1][this.ij[1] - 1] > this.zMaxMin) {
		    this.elle = this.idir % 2;
		    this.ij[this.elle] = this.sign(this.ij[this.elle], this.l1[this.k - 1]);
		    continue;
		}
		this.elle = 0;
		return false;
	    }
	    if (this.idir != this.nxidir) {
		this.nxidir++;
		this.ij[this.elle] = this.l1[this.k - 1];
		this.k = this.nxidir;
		this.elle = 1 - this.elle;
		this.ij[this.elle] = this.l2[this.k - 1];
		if (this.nxidir > 3) {
		    this.nxidir = 0;
		}
		continue;
	    }

	    if (this.ibkey != 0) {
		return true;
	    }
	    this.ibkey = 1;
	    this.ij[0] = this.icur;
	    this.ij[1] = this.jcur;
	    if (this.routineLabel020()) {
		continue;
	    }
	    return false;
	}
    }

    // -------------------------------------------------------
    // "routineLabel200" corresponds to a block of code
    // starting at label 200 in Synder's subroutine "GCONTR".
    // It has return values 0, 1 or 2.
    // -------------------------------------------------------
    private int routineLabel200(final boolean workSpace[]) {
	while (true) {
	    this.xy[this.elle] = (1.0 * this.ij[this.elle]) + this.intersect[this.iedge - 1];
	    this.xy[1 - this.elle] = 1.0 * this.ij[1 - this.elle];
	    workSpace[(2 * (((this.xSteps * (((this.ySteps * this.cntrIndex) + this.ij[1]) - 1)) + this.ij[0]) - 1)) + this.elle] = true;
	    this.drawKernel();
	    if (this.iflag >= 4) {
		this.icur = this.ij[0];
		this.jcur = this.ij[1];
		return 1;
	    }
	    this.continueContour();
	    if (!workSpace[(2 * (((this.xSteps * (((this.ySteps * this.cntrIndex) + this.ij[1]) - 1)) + this.ij[0]) - 1)) + this.elle]) {
		return 2;
	    }
	    this.iflag = 5; // 5. Finish a closed contour
	    this.iedge = this.ks + 2;
	    if (this.iedge > 4) {
		this.iedge = this.iedge - 4;
	    }
	    this.intersect[this.iedge - 1] = this.intersect[this.ks - 1];
	}
    }

    // -------------------------------------------------------
    // "crossedByContour" is true iff the current segment in
    // the grid is crossed by one of the contour values and
    // has not already been processed for that value.
    // -------------------------------------------------------
    private boolean crossedByContour(final boolean workSpace[]) {
	this.ii = this.ij[0] + this.i1[this.elle];
	this.jj = this.ij[1] + this.i1[1 - this.elle];
	this.z1 = this.z[this.ij[0] - 1][this.ij[1] - 1];
	this.z2 = this.z[this.ii - 1][this.jj - 1];
	for (this.cntrIndex = 0; this.cntrIndex < this.ncv; this.cntrIndex++) {
	    final int i = (2 * (((this.xSteps * (((this.ySteps * this.cntrIndex) + this.ij[1]) - 1)) + this.ij[0]) - 1)) + this.elle;

	    if (!workSpace[i]) {
		final double x = this.cv[this.cntrIndex];
		if ((x > Math.min(this.z1, this.z2)) && (x <= Math.max(this.z1, this.z2))) {
		    workSpace[i] = true;
		    return true;
		}
	    }
	}
	return false;
    }

    // -------------------------------------------------------
    // "continueContour" continues tracing a contour. Edges
    // are numbered clockwise, the bottom edge being # 1.
    // -------------------------------------------------------
    private void continueContour() {
	int local_k;

	this.ni = 1;
	if (this.iedge >= 3) {
	    this.ij[0] = this.ij[0] - this.i3[this.iedge - 1];
	    this.ij[1] = this.ij[1] - this.i3[this.iedge + 1];
	}
	for (local_k = 1; local_k < 5; local_k++) {
	    if (local_k != this.iedge) {
		this.ii = this.ij[0] + this.i3[local_k - 1];
		this.jj = this.ij[1] + this.i3[local_k];
		this.z1 = this.z[this.ii - 1][this.jj - 1];
		this.ii = this.ij[0] + this.i3[local_k];
		this.jj = this.ij[1] + this.i3[local_k + 1];
		this.z2 = this.z[this.ii - 1][this.jj - 1];
		if (((this.cval > Math.min(this.z1, this.z2)) && (this.cval <= Math.max(this.z1, this.z2)))) {
		    if ((local_k == 1) || (local_k == 4)) {
			final double zz = this.z2;

			this.z2 = this.z1;
			this.z1 = zz;
		    }
		    this.intersect[local_k - 1] = (this.cval - this.z1) / (this.z2 - this.z1);
		    this.ni++;
		    this.ks = local_k;
		}
	    }
	}
	if (this.ni != 2) {
	    // -------------------------------------------------
	    // The contour crosses all 4 edges of cell being
	    // examined. Choose lines top-to-left & bottom-to-
	    // right if interpolation point on top edge is
	    // less than interpolation point on bottom edge.
	    // Otherwise, choose the other pair. This method
	    // produces the same results if axes are reversed.
	    // The contour may close at any edge, but must not
	    // cross itself inside any cell.
	    // -------------------------------------------------
	    this.ks = 5 - this.iedge;
	    if (this.intersect[2] >= this.intersect[0]) {
		this.ks = 3 - this.iedge;
		if (this.ks <= 0) {
		    this.ks = this.ks + 4;
		}
	    }
	}
	// ----------------------------------------------------
	// Determine whether the contour will close or run
	// into a boundary at edge ks of the current cell.
	// ----------------------------------------------------
	this.elle = this.ks - 1;
	this.iflag = 1; // 1. Continue a contour
	this.jump = true;
	if (this.ks >= 3) {
	    this.ij[0] = this.ij[0] + this.i3[this.ks - 1];
	    this.ij[1] = this.ij[1] + this.i3[this.ks + 1];
	    this.elle = this.ks - 3;
	}
    }

    // -------------------------------------------------------
    // "contourPlotKernel" is the guts of this class and
    // corresponds to Synder's subroutine "GCONTR".
    // -------------------------------------------------------
    private void contourPlotKernel(final boolean workSpace[]) {
	int val_label_200;

	this.l1[0] = this.xSteps;
	this.l1[1] = this.ySteps;
	this.l1[2] = -1;
	this.l1[3] = -1;
	this.i1[0] = 1;
	this.i1[1] = 0;
	this.i2[0] = 1;
	this.i2[1] = -1;
	this.i3[0] = 1;
	this.i3[1] = 0;
	this.i3[2] = 0;
	this.i3[3] = 1;
	this.i3[4] = 1;
	this.i3[5] = 0;
	this.prevXY[0] = 0.0;
	this.prevXY[1] = 0.0;
	this.xy[0] = 1.0;
	this.xy[1] = 1.0;
	this.cntrIndex = 0;
	this.prevIndex = -1;
	this.iflag = 6;
	this.drawKernel();
	this.icur = Math.max(1, Math.min((int) Math.floor(this.xy[0]), this.xSteps));
	this.jcur = Math.max(1, Math.min((int) Math.floor(this.xy[1]), this.ySteps));
	this.ibkey = 0;
	this.ij[0] = this.icur;
	this.ij[1] = this.jcur;
	if (this.routineLabel020() && this.routineLabel150()) {
	    return;
	}
	if (this.routineLabel050()) {
	    return;
	}
	while (true) {
	    this.detectBoundary();
	    if (this.jump) {
		if (this.ix != 0) {
		    this.iflag = 4; // Finish contour at boundary
		}
		this.iedge = this.ks + 2;
		if (this.iedge > 4) {
		    this.iedge = this.iedge - 4;
		}
		this.intersect[this.iedge - 1] = this.intersect[this.ks - 1];
		val_label_200 = this.routineLabel200(workSpace);
		if (val_label_200 == 1) {
		    if (this.routineLabel020() && this.routineLabel150()) {
			return;
		    }
		    if (this.routineLabel050()) {
			return;
		    }
		    continue;
		}
		if (val_label_200 == 2) {
		    continue;
		}
		return;
	    }
	    if ((this.ix != 3) && ((this.ix + this.ibkey) != 0) && this.crossedByContour(workSpace)) {
		//
		// An acceptable line segment has been found.
		// Follow contour until it hits a
		// boundary or closes.
		//
		this.iedge = this.elle + 1;
		this.cval = this.cv[this.cntrIndex];
		if (this.ix != 1) {
		    this.iedge = this.iedge + 2;
		}
		this.iflag = 2 + this.ibkey;
		this.intersect[this.iedge - 1] = (this.cval - this.z1) / (this.z2 - this.z1);
		val_label_200 = this.routineLabel200(workSpace);
		if (val_label_200 == 1) {
		    if (this.routineLabel020() && this.routineLabel150()) {
			return;
		    }
		    if (this.routineLabel050()) {
			return;
		    }
		    continue;
		}
		if (val_label_200 == 2) {
		    continue;
		}
		return;
	    }
	    if (++this.elle > 1) {
		this.elle = this.idir % 2;
		this.ij[this.elle] = this.sign(this.ij[this.elle], this.l1[this.k - 1]);
		if (this.routineLabel150()) {
		    return;
		}
	    }
	    if (this.routineLabel050()) {
		return;
	    }
	}
    }

    /*
     * attach the segment to the existing one or create a new segment
     * 
     * iflag == 1 means Continue a contour iflag == 2 means Start a contour at a
     * boundary iflag == 3 means Start a contour not at a boundary iflag == 4
     * means Finish contour at a boundary iflag == 5 means Finish closed contour
     * not at boundary iflag == 6 means Set pen position
     */
    private void contourSegment(final double prevU, final double prevV, final double u, final double v, final double z) {
	int i = 0;
	for (i = 0; i < this.contours.size(); i++) {
	    final Contour c = this.contours.get(i);
	    if (!c.finished) {
		final Point p = c.data.get(c.data.size() - 1);
		if ((p.getX() == prevU) && (p.getY() == prevV)) {
		    c.data.add(new Point(u, v));
		    if (this.iflag == 4) {
			c.finished = true;
		    } else if (this.iflag == 5) {
			c.finished = true;
			c.closed = true;
		    }
		    break;
		}
	    }
	}
	if (i == this.contours.size()) { // start a new contour
	    final Contour c = new Contour();
	    c.z = z;
	    c.data.add(new Point(prevU, prevV));
	    c.data.add(new Point(u, v));
	    this.contours.add(c);
	}
    }
}
